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We present a technique to resolve a Gaussian density matrix and its time evolution through known 
expectation values in position and momentum. Further we find the full spectrum of this density 
matrix and apply the technique to a chain of harmonic oscillators to find agreement with conformal 
field theory in this domain. We also observe that a nonconformal state has a divergent entanglement 
entropy. 
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I. INTRODUCTION 

Entanglement is today considered a fundamental re- 
source in nature when it comes to quantum computation 
and information [[J, and measures of entanglement has 
become a major field of research. In particular, entan- 
glement in condensed-matter systems and the entangle- 
ment's critical behavior is well investigated @,0,HU0|- 
In these terms, the entanglement entropy is an analyt- 
ically well suited tool for investigating the properties 
of ground states in condensed-matter systems. In this 
paper, we focus on bososnic states with Gaussian wave 
functions, and the entanglement properties of the ground 
state of a simple harmonic chain, which belongs to this 
class 0. 

We consider the notion of entanglement entropy, that 
is, considering a quantum system denoted C in a pure 
state with wave function \ip), we trace out some portion 
B to obtain the density matrix of the remaining space A 
as pa = Tr g | ip) (if) | . Then the entanglement of A with 
respect to B is well defined by the entropy (measured in 
ebits) of the reduced density matrix 0, Q ; 



remaining space with respect to the rest is 



Sa = -Trp^log 2 p A . 



(1) 



The measure of entropy in units of ebits is customary and 
we will use logarithms base two throughout the paper. 
This procedure is well established, and works well for all 
cases where the entire strip C is in a pure state, though 
entanglement measures for mixed states are still incom- 
plete. Most work has been focused on this entanglement 
in spin chains, but we will focus on a one-dimensional 
bosonic strip |J El ■ 

At critical points in a parameter space we have scale 
and translational invariance, and thus expect the theory 
to be conformally invariant ^lj, and one can use this 
fact to efficiently detect critical systems J5§. As was com- 
puted by Holzhey, Larsen, and Wilczek |12j . conformally 
invariant systems in 1+1 dimension can be considered as 
a string of length A of which we trace out some fraction 
1 — <t € [0,1], and then the entanglement entropy of the 



S(a) = 



G 



dog 



Sm(7T<7j 

7T6 



(2) 



Here c and c are the holomorphic and anti holomorphic 
charges respectively, e is some cutoff parameter that 
we will consider arbitrary. When considering the limit 
a <C 1, the formula reduces simply to S ~ log oA, which 
has been a matter of keen interest 0- However, we will 
focus on any er, in particular, when keeping the A con- 
stant (J2J) provides a very specific signature of a confor- 
mally invariant system. Thus this formula presents two 
independent (as long as e is considered arbitrary) signa- 
tures of a finite conformal system. First the logarithmic 
divergence of the entropy as a is held constant while A 
increases, and second the characteristic log sin signature 
when A is constant and a varies. 

Another, more trivial, measure of the entanglement of 
a reduced density matrix is the product state identifica- 
tion 



M 



E M = 1-Trp^ 



M > 2, 



(3) 



which is zero for a product state, and unity for a maxi- 
mally entangled state. This measure is equivalent to the 
Renyi entropy |l3j and is not well suited for much more 
than to single out a pure state, as with increasing M any 
entangled state will converge to zero in this measure, 



lim Em 



1 entangled state 
product state. 
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The density matrix contains all information of any sys- 
tem in a mixed or pure state, and computing vital physi- 
cal information on any such system is mostly determined 
by the eigenvalues of the density matrix. Therefore there 
is a great need to compute these in an efficient way. In 
particular, considering a bosonic system in one spatial 
dimension in its ground state, it can be modeled as a 
harmonic chain. The quantum correlations as measured 
by the entanglement are nonzero, as will be shown. Thus 
the vacuum itself is entangled, which is a highly intrigu- 
ing result. Our main purpose in this paper is twofold, 
first we wish to demonstrate how to compute this entan- 
glement from the vacuum ground state, and second we 
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wish to see how and if the conformal invariance identi- 
fied by Eq. J2J arises in this case. The free boson is 
known to have a central charge c = c = 1, so the pref- 
actor in S(a) reduces to simply 1/3. It turns out that 
we indeed have a conformally invariant ground state, as 
identified by the two conformal signatures provided. As 
we will demonstrate in Sec. IIVI the entropy diverges in 
the massless limit, but it nevertheless seems conclusive 
that the theory is conformally invariant below a certain 
threshold mass. 

In this paper we will in Sec. HTIcompute how a Gaussian 
density matrix can be recovered from expectation values 
in position and momentum, and how this can be used to 
compute the time evolution of the density matrix. This 
can also be computed in terms of Weyl representation, 
though we apply a more explicit representation here. In 
Sec. IHII we compute the two entanglement measures writ- 
ten in this section from the Gaussian density matrix, and 
finally in Sec. IIVI we apply the formulas to the ground 
state of a harmonic chain. 



Q and P are symmetric, while all are real. Furthermore, 
since the system in translationally invariant we can as- 
sume Qij = f(\i — j\) for some function /, and simi- 
larly with the other matrices. Thus Q, P, and S are 
Toeplitz matrices, potentially simplifying numerical com- 
putations. Toeplitz matrices are known to be central also 
in the study of quantum spin chains These three 

matrices can be assumed known in a given model of a 
bosonic quantum system. We refer to the set of variables 
in the density matrix as 9 = {A, C, d}, and the expecta- 
tion value matrices as S = {Q, P, S, (q), (p)}. A simple 
count shows that the two sets both have 2D 2 + 3D de- 
grees of freedom, and thus the two sets may contain equal 
amounts of information. The expectation value matri- 
ces can be computed for a given density matrix through 
a straightforward calculation that computes expectation 
values of an operator O as 



(6) =Tr6p(q,q') 



dqdq'S(q-q')6p(q,q'). 



II. GAUSSIAN DENSITY MATRICES 

A. Recovery of the density matrix from 
expectation values 

It is known that a Gaussian state in a quantum har- 
monic oscillator potential will evolve in such a way that 
the Gaussian shape is preserved at all times. However, 
the different parameters of the Gaussian state may also 
evolve in time so that the exact appearance of the density 
matrix may be difficult to predict. Consider a Gaussian 
density matrix of N particles with positions qi, 



/ det (A'-C) 

p{q, q) = \l rn ex P 



x exp 
+ d*q', 



(A'ij — C-j) d'j 



L y ('/<•! ..,'/.. + '/!• !, ',</',) + '...'/', + <kqi 

(4) 



summing over repeated indices. Here A is a positive, 
symmetric N x N matrix, while C is a Hcrmitian N x 
N matrix, and d is an A^-dimensional vector. We use 
A', C, and d! to denote the real parts of A, C, and 
d, respectively. We further denote the imaginary parts 
such that A = A' + iA", etc. Positions and momenta are 
real valued. The matrix A' — C must be invertible and 
positive in order to have an normalizable density matrix. 
Also, in order to keep p positive, that is that (ip\p\ip) > 
for any one must have A' > 0. Hence follows that 
even C must be positive. Next we define three matrices 
of variances in position and momentum, 



Pi, 



(QiOj) - 
(PiPj) 



(5) 



s%j = ^(qiPj +Pjq*) - (ft>(pj>- 



This means that the expectation value matrices S[0] can 
be computed in terms of the parameters in the density 
matrix, 



Q = \(A'-C')-\ 

P = A - (A - C)Q(A — C) T , 
S = -Q{A" + C"), 
(q) = 2Qd', 

(p) = -2{A" - C")Qd' + d". 



(6) 



Note that since A' — C is invertible, Q is well defined 
invertible. This set of matrix equations can then be in- 
verted to yield 9 [3], 

A' = P+^Q- 1 - S T Q- 1 S, 
A" = - 1 -(S T Q- 1 + Q- 1 S), 



C = P- -Q- 1 - S T Q- 1 S, 



(7) 



d'^-Q-^q), 
d" = {A" + C"){q) + {p). 

We note that S = implies both A and C to be real, 
and that in a state where (q) = (p) — 0, d = 0. This 
way it is simple to recover the density matrix given S. 
A strategy to find the nonlinear time evolution of the 
density matrix is sketched in Fig. ^ one would use the 
fact that the time evolution of S is much simpler to find 
than that of 9. And since the formulas above enables us 
to switch between the two sets in an easy manner, we can 
take the time evolution via 3 instead of taking it directly 
on 9. 
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9(0) _ 
Eq.© 



6(f) 



Eq.0 



"/q\ Linear evolution . 



FIG. 1: Outline of our strategy to find the time evolution of 
Q(t) by going through E(i). 



B. Time evolution 

The Gaussian shape of the density matrix is preserved 
under time evolution governed by the Lagrangian 

£ = ^ iinin — Qn^nmQm + £,nQn) , (8) 

with sum over repeated indices. The conjugate momen- 
tum is p n = q n , while O is symmetric and £ is a real 
three- vector describing external forces. The equations of 
motion become q n + 0, nm q m = £„, which implies four 
coupled differential equations for the time evolution of 
the expectation value matrices 



Q = T(S), 
P = -T(flS), 



S 



where we have defined the symmetrizing operator 
T{A) = A + A T . These relations specify the time evo- 
lution of the system, given some initial condition. Also, 
combining the first and third equation gives 



Q + 2Q--{Q,tl} 



0, 



or, in the case of a diagonal fi, fij,- = SijUi, this reads 



Qn + 2Q 



= 0. 



III. ENTANGLEMENT MEASURES FOR 
GAUSSIAN DENSITY MATRICES 



with the mentioned requirement on the density matrix 
that A' — C is positive. We will omit the two latter 
parameters in the argument list of p° when there is no 
risk of confusion. The eigenvalue equation of this density 
matrix is 

J dx'p°(x,x')^ n (x') = A„* n (x). 



Also, the density matrix obeys 



p°(x + x , x' + .t ) = e id " x p°(x, x'; rj, 0)e 



— id' x 



(ii) 



where xq — d'/l — rj. This means that the scaled eigen- 
function ^ n (x) = e ld ( x ~ x ° > '$> n {x — xq) also is an eigen- 
vector of p (x, x') with the same eigenvalue as 4 , n (x), and 
hence the eigenvalues are independent of the displace- 
ment d. Furthermore, Eq. shows that any traces 
Tr p M are invariant under d. 

Now, to find the eigenvalues of the single-particle den- 
sity matrix when d = 0, consider the Green's function 
Q(z,z';t) = (z\e~ rH \z') , with H being the single har- 
monic oscillator Hamiltonian H = —p 2 /2 + \uz 2 , and 
the states \z) the position eigenstates of the quantum 
harmonic oscillator. The eigenvalue result for this matrix 
is, with {z\^) n = and \^>} n , the energy eigenstates 



dz' g(z, z'; r)*„(/) = e-^ n+1 ' 2 ^ n {z). 



(9) Furthermore, Q must solve the initial value problem 



-^-Hjg(z,z';r)=0, 
lim Q{z, z'\ t) = S(z — z'), 



which, through a scaling of variables, x — Zy/ui coth(wr), 
has the solution 



Q(x,x';t) 



2tt sinh(ti>r) 



1 



x exp 



{x 2 + x' 2 ) + 



cosh(o;r) 



(12) 



There also exists a complex solution, which we disregard 
as unphysical. This solution (|12J) means that if we iden- 
tify X] — 1/ cosher), we can conclude that p° has the 
infinite set of eigenvalues 



The eigenvalues of Eq. 10} can be found explicitly as 
we will show in this section. To this end, consider what 
we will refer to as the single-particle density matrix 



p (x, x') r],d) = e 



- p -d' 2 /l-V 



1 



V 



x exp 



— ( x 
2 



(x 2 + x' 2 ) + r]xx' + dx + d*x' 



(10) 



Here d = d! + id" is a complex number, while r\ is real 
in the range [0,1). The latter constraint is consistent 



A„ = ^ujcoth(LUT)e-^ T{n+1/2) = A C (13) 

ne {0,1,2,...}, 



where £ = e~ UT = T]/(l + y 1 — rf ) . 

Hence we have identified the eigenvalues and eigen- 
functions of our density matrix in the case d = 0, which 
is sufficient to know the eigenvalues also in the rf ^ 
case. 

Now, let us again turn to our general density matrix 
(0} , and show that this can be transformed into the sim- 
ple form of Eq. I|10|) 14]. In general, the matrix A" 
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contributes nothing to traces of the density matrix, and 
A can be considered real. Take an orthogonal matrix O 
that diagonalizes a = O t A'O and create the matrix y/a 
which is diagonal and takes the square roots of A's eigen- 
values on its diagonal. A is positive by assumption, so 
this is well defined and real. Then we find the orthogonal 
matrix O that diagonalizes the matrix f\ defined as 



V 



(14) 



Note that both O and O are real, since C is hermitian and 
A is symmetric. The coordinates and vectors d transform 
as 

q = OV^O T q, 
d = dVF T T d, 



a result that agrees exactly with that found by Srednicki 
IS' 

We can proceed to find exact formulas for the entan- 
glement measures Em = 1 — Tr p M where M > 2, for 
which we find that 



N fl -f ) M 

t=l 



M 



(18) 



These two entanglement measures are now quite 
straightforward to compute. Since Q is a Toeplitz ma- 
trix, the inversion involved can be done efficiently with 
existing linear algebra packages. Also the computation 
involves two diagonalizations of real, symmetric N x N 
matrices, which is also numerically straightforward and 
efficient. We will focus primarily on the entropy due to 
its analytical usability. 



which means that the density matrix J3J may be rewrit- 
ten as 



IV. APPLICATION 



N 



(15) 



i=l 



rji being the diagonal elements of fj and p(q, q') = 
p(q,q'). In other words, the density matrix Q| may be 
expressed as a product of single-particle density matrices 
ifTHjl . Thus each of the product terms in Eq. lfT5*|> gives 
the infinite set of eigenvalues l|13l) , and we may write the 
eigenvalues for the total matrix as 



N 



N 



(16) 



Here we have defined Ao 



Here we have let 

Ao i— > \q' , Tf i — > Tji) and £ h- > £j from Eq. (|13|) where the 
index i refers to particle number. Ao can be calculated 
from the normalization condition 



n D = l 



giving A = — This sets the stage for the 

calculation of entropy and traces over the density matrix. 



Entanglement measures 

Having obtained the density matrix eigenvalues, the 
entropy S = — Tr p\ogp is easily calculated, 



6 log 6 



ni — l njy — 1 

N r 

log (!-&) + 



E 



i-6 



(17) 



We now apply the above formalism to the quantum 
Klein-Gordon field </>„, defined by the Lagrangian © 
with £ = and f2y = J^k on iV lattice points with 
periodic boundary conditions [l6j . The lattice constant 
is denoted a, and the system size is thus A = aN. This 
field has the Fourier expansion in bosonic creation and 
annihilation operators afc and a\. respectively; 



T — ( 

4^ 2Au k V 



and conjugate field 




fl t gi(fcn-w fc t) 



(20) 



In this discrete field theory the dispersion relation is 



"2 



— sin 2 (fc/2) + K 2 



and the sum is over all allowed wave vectors k in this 
space. We have the nonzero commutation rules (at equal 
time) 

[4>(x),tt(x')] — iS(x — x'), 

[Ofe, <4j = S k .k': 

indicating the harmonic oscillator nature of the system, 
with positions <J) n and momenta n n . When rescaling this 
theory it is imperative to keep A constant in order to 
maintain any scaling invariance and hence any conformal 
invariance in the model. The S-matrices of the ground 
state are the Toeplitz matrices 



(m—n) 



Pn 

Sri 



2N ^ 

k 

o, 



Lu k e 



ik(m- if ) 



(21) 
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FIG. 2: E2 as function of k for three different system sizes. 
The transition from an entangled state in the massless limit 
to a product state in the massive limit is clearly shown. Here, 
and in Fig. |3] n is in units of inverse length 



FIG. 3: Smax with n for the same system sizes as in Fig. |5] 
The three lower lines are a = 1, while the divergent lines are 
a — 0. Note that the two lines with different a and same N 
diverge from each other at roughly the same n regardless of 
the system size N. 



with the sums running over all lattice points numbered 
i, and k = 2iri/N . Q and P are thus essentially finite 
Fourier transforms of Wj. 1 . Now, since S is vanishing, 
we can conclude that both A and C are real. As shown 
this is all we need to recover the density matrix, and 
calculating the traces. 

In the limit of a massless theory, k — > 0, Q n will diverge 
due to the k = (zero mode) term in the sum jZH - We 
can exclude this zero mode from the sum by summing 
over lattice point m using k = 7r(2m + a) /N and choose 
a = (1) to compute with (without) the zero mode. This 
corresponds to using (anti)periodic boundary conditions 
on the field 4>. 

When finding the entropy we define some of the points 
as "inside" (region A) and some as "outside" (region B) , 
and trace over the outside to find the geometric entropy 
of the inside region. Formally, this amounts to calculat- 
ing the 3 matrices for the crN inside oscillators, a be- 
ing the fraction of the entire system that A constitutes. 
When the Ss are known, we can compute the density ma- 
trix and the entropy for the state. Clearly, the entropy 
is symmetric with respect to interchange of A and B, 
so any entanglement measure will be symmetric around 
a = 1/2. Hence we expect a maximal entanglement at 
this half size, and we can consider the half size entangle- 
ment with respect to k and N. Results for E2 at half 
size is shown in Fig. [21 Hence we conclude that the mas- 
sive k — ► 00 system in nonentangled, while we have a 
transition to a maximally entangled case in the massless 
system k — > 0. Also, the entanglement is larger for larger 
systems, and the transition occurs at a larger mass in 
a larger system. The correlations are greater in a large 
system, and the inertia of the mass k must be larger to 
prevent them. 

We now consider the entanglement entropy. The half 
size measure for both a = and a = 1 are shown in 




N 



FIG. 4: S m ax with iV for various k. The straight lines shown 
all rise with central charge c = 1/3. The bisected points are 
the a = points, from top to bottom, k = 10~ 4 , k = 10~ 2 , 
K = 1, and k = 10 2 . The circles are for a = 1 and k = 0. 



Fig. |31 We see that when the zero mode is included, 
the entropy diverges in the massless limit, while without 
the zero mode it converges to some system dependent 
value. In Fig. 0] we see how the half size entropy di- 
verges logarithmically with system size for small k. Most 
remarkably, this divergence occurs regardless of the zero 
mode, and the correct scaling factor c/3 is reproduced. 

Finally, we look at the second feature of the conformal 
signature, namely the log sin shape of a finite system. For 
the a = case this is shown in Fig. and we again see 
a good characteristic of the conformal system in a mass- 
less system. For the massive system, the entanglement 
saturates, and fits the signature only at small er, which is 
also observed earlier in noncritical quantum spin chains 
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FIG. 5: S(a) for some values of n with a — 0. Here, TV = 
100(+) and TV = 10(O)- From top to bottom, k = 10" 4 , 
K — 10 -2 , K — 1, and k = 10 2 . The full lines are /(<x) = 
| log 2 (sin 7T(j) + a, with a chosen to fit the lines at the ends. 
For large ac, the entanglement saturates and does not obey 
the conformal signature. In the massless limit the conformal 
signature is obeyed, even at TV = 10. 
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FIG. 6: The largest entropy contributions in the sum S = 
~}2 n s n, with ordered terms. The graph shows the terms for 
system size TV = 100, a = (+), and a = 1 (x) and also 
system size TV = 25, a — (*), and a = 1 (□). All values are 
for an essentially massless theory with k = 10 -3 . 

0. When the zero mode is omitted, however, the log sin 



signature is not present, as the system is identical to the 
a = case for some massive k, a state that is not confor- 
mally invariant. But then it is nevertheless notable that 
the logarithmic divergence with system size still fits the 
conformal theory although the fixed-size signature does 
not, and the state is not conformally invariant. 

It is valuable to note that the conformal signature in 
the entanglement entropy is present even in small sys- 
tems, such as TV = 10. This has been seen in quantum 
spin chains earlier, there enabling efficient identification 
of criticality . 

As a final feature, we investigate which modes, or 
terms, contribute to the entanglement entropy in the ex- 
pansion S = J2 n s n where the terms arc defined by the 
sum (11711 . The individual terms are shown in Fig. [SI 
where we see that for the conformally invariant case, the 
eigenvalues fall of faster than exponentially, which indi- 
cates that only very few terms in the sum are needed to 
compute the entropy. Indeed, to compute the entropy to 
within an error ±10~ 6 one needs only eight out of a pos- 
sible 50 terms in the expansion. Also, the figure shows 
that for the nonconformal state without the zero mode, 
the eigenvalues are paired since in a not scale invariant, 
though translationally invariant state, any mode with a 
nonzero impulse will be degenerate with another state of 
opposite impulse. 



V. CONCLUSION 

We have shown how to compute the entanglement en- 
tropy of a chain of bosonic harmonic oscillators, using 
scalings and rotations of the density matrix to put it in a 
single-particle form. The results show what we call con- 
formal signatures in the massless limit. Moreover, the 
results show that even a nonconformal state can show 
a logarithmic divergence as predicted by conformal field 
theory, but not the log sin signature that the author be- 
lieves to be a uniquely conformal feature. 
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